  #include <stdio.h>
  #include <stdlib.h>
  #include <math.h>

  #include <common.h>
  #include <parse_config_file.h>

  double r_rand( double l, double h ) {
  
    return l + drand48()*(h-l);
  }
  
  void make_plummer_model( SIM_DATA *s ) {
    
    int i;
    double M,r,t,p,x,y,z,v,vx,vy,vz,xx,yy;
    long N = s->N;
    for( i=0; i< N; i++ )
    {
      M = 1./N;
      r = 1.0 / sqrt( pow(r_rand(0,1), -2.0/3.0) - 1.0);
      t = acos( r_rand(-1,1) );
      p = r_rand(0,2*M_PI);
      
      x = r * sin(t) * cos(p);
      y = r * sin(t) * sin(p);
      z = r * cos(t) ;
      
      xx=0,yy=0.1;
      while( yy > xx*xx*pow(1-xx*xx,3.5) )
      {
        xx = r_rand(0,1);
        yy = r_rand(0,0.1);
      }
      
      v = xx * sqrt(2.0) * pow( 1.0 + r*r,-0.25);
      t = acos( r_rand(-1,1) );
      p = r_rand(0,2*M_PI);
      
      vx = v * sin(t) * cos(p);
      vy = v * sin(t) * sin(p);
      vz = v * cos(t) ;
      
      s->X[4*i+0] = x; s->X[4*i+1] = y; s->X[4*i+2] = z; s->X[4*i+3] = M;
      s->XDOT[4*i+0] = vx; s->XDOT[4*i+1] = vy; s->XDOT[4*i+2] = vz;
    }
    
  }
 
  int main( int argc, char *argv[] ) {
	  
	if( argc != 3 ) {
      fprintf( stderr, "Usage: ./makeplummer <config file> <output file\n " );
      return 1;
	}
		
    SIM_DATA sim0;  

    FILE *in_config;
    in_config = fopen( argv[1], "r" );	
	
	//Careful with this, it will overwrite!
	FILE *outfile;
	outfile = fopen( argv[2], "w" );
	
  // Load the configuration file of the simulation
    parse_config_file( in_config, &sim0 );
  
  // Allocate neccesary memory
    allocate_arrays( &sim0 );
  
  // Check everything is in order
    check_simulation_params( &sim0 );
	
  // Build the plummer model
	make_plummer_model( &sim0 );
	
  // Dump the output
	dump_phasecoordinates( outfile, &sim0 );
	
  // Close the streams
	fclose( in_config );
	fclose( outfile   );
	
	return 0;
  }
